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Abstract 

We point out a new mechanism which can lead to mean field type be- 
haviour in nonequilibrium critical phenomena. We demonstrate this mech- 
anism on a two-dimensional model which can be understood as a stochas- 
tic and non-conservative version of the abelian sandpile model of Bak et 
al. [Qj. This model has a second order phase transition whose critical be- 
haviour seems at least partly described by the mean field approximation for 
percolation, in spite of the low dimension and the fact that all interactions 
are of short range. Furthermore, the approximation obtained by replacing 
the lattice by a Bethe tree is very precise in the entire range of the control 
parameter. 
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Although mean field theory is one of the most useful and most commonly ap- 
plied approximations for spatially extended complex systems, it has an important 
drawback: near critical points it is usually not only quantitatively but also quali- 
tatively wrong, in the sense that it predicts wrong critical exponents and scaling 
functions. There are only few known exceptions to this. They involve systems 
with infinitely long range interactions, systems in high dimensions, and systems 
defined on "lattices" without loops ("Bethe tree" lattices). 

The reason why mean field theory gives wrong results is easily understood quali- 
tatively. By neglecting correlations, it also neglects the feedback of the local order 
parameter value s(x) onto itself. Since a system at its critical point is very sen- 
sitive to any influence, these feedbacks are enhanced and modify its behaviour in 
an essential way. The above exceptions avoid this problem by either disallowing 
feedback loops (Bethe lattices) or by suppressing the effect of short loops in favour 
of long loops which cover essentially the whole lattice. 

In this note we want to point out a novel mechanism for generating mean field 
behaviour. This mechanism can only be effective in non-equilibrium systems, but 
we conjecture that it should occur there rather commonly. In equilibrium systems, 
the above-mentioned feedback is always positive due to detailed balance. In non- 
equilibrium systems, however, the feedback can also be negative. If the system 
reacts to the feedback in a strongly non-linear way, the effects of individual loops 
will thus not be additive, but might cancel. In general, this cancellation will 
of course not be complete, and the net effect of the feedback will be similar to 
random noise. Otherwise said, the order parameter s(x) will be "confused" by the 
influences of the different feedback loops. This confusion will lead to a fast loss 
of memory which in turn implies that the mean field approximation can become 
correct. 

Specifically, we shall illustrate this with a model which is inspired by the sandpile 
model []]]. It is also abelian in the sense of @, but the evolution of avalanches 
is stochastic and non- conservative. As a consequence, its critical behaviour is not 
self- organized, but it has a critical point in the ordinary (co-dimension 1) sense. 
A model which gives identical spatial structures was first introduced by one of us 
in || (model nr. 5 in the appendix of that paper). Similar (but somewhat more 
complicated) models were also studied in Q. Although we cannot prove analyti- 
cally that the critical behaviour of our model is of mean field type, our numerical 
results suggest very strongly that at least two of the three independent critical 
exponents have their mean field value, and that scaling functions are identical. 

Our model is defined on a two-dimensional square lattice (generalizations to higher 
dimensions are obvious), and time is discrete. At each lattice site we have a "spin" 
Zi, which can take any non-negative integer value, but only the values Zi = and 
Zi — 1 are "stable". If Zi becomes > 1 during the evolution, it "topples" at the 



2 



next time step. The toppling rule is 



i-2 



(1) 



for the site which topples, and 




Zj + 1 

Z 3 



with probability p 
with probability 1 — p 



(2) 



for each of its four neighbours. Notice that each neighbour j has the same chance 
p to get its spin increased, independently of what happens at the other neighbours. 
Thus the sum J2i z % fluctuates (during each toppling, it can change by any value 
between -2 and +2), but in the average each toppling event causes J2i %i to decrease 
by 4p — 2. The critical point is exactly where this vanishes, p c = 1/2. For later 
use we define e = 1/2 — p, and g = (z). 

As in the sandpile model, the dynamics consists of two rules: 

(i) If all sites are stable, a site i is chosen randomly, and Zi as well as time t are 
increased by one unit. In the following, we call this an "event". 

(ii) If at least one site is unstable, the above toppling rule is applied simultane- 
ously at all unstable sites. If some Zi > 4 so that it has to topple repeatedly to 
become stable, then all these topplings are also done simultaneously. After this, t 
is increased by 1. This is repeated until all sites are stable again, after which rule 
(i) is applied again. 

Our most extensive studies were done for p < p c , where all avalanches are finite 
(with probability one), even on an infinite lattice. In this case we can use periodic 
boundary conditions, which greatly reduces finite-size effects. Some simulations 
were done also at p = p c . There g would not decrease during an avalanche, leading 
to problems with infinite avalanches when periodic boundary conditions are used 
[[ In this case, we used the same kind of open b.c. as in the sandpile model: if a 
neighbour of site i is outside the boundary, we just disregard it in eq.(|2|). 

In a mean field theory for this model, we assume that a toppling site "forgets" 
after the next time step that it has toppled, and Zj is replaced by an identically and 
independently distributed random variable which takes the value 1 with probability 
g, and with probability 1 — g. This implies in particular that there are no 
correlations, c(i — j) = (ziZj) — g 2 = 0. Correlations measured in simulations are 
shown in fig. 1. We see that they are indeed small and, what is more important, 
they decay very rapidly with distance. A similarly rapid decay was observed also 
in the sandpile model || . In that model it was indeed proven that the correlations 
decay as ||i — j\\~ 4 ||. 

Notice that we assume that this loss of memory due to confusion by the topplings 
of the other neighbours happens only after the next time step. So each toppling site 

1 Serious problems arise only when p > p Cl but we expect very slow behaviour and numerical 
instabilities already at the critical point. 
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Figure 1: Correlations measured in simulations, plotted against the square of the Eu- 
clidean distance. For distances > 4, all correlations were compatible with zero. 



(except the very first of the avalanche which we call the root) has three neighbours 
which it can cause to become unstable with probability a = gp, while the fourth 
neighbour (its father) will be induced to topple with a different probability a'. 
The simplest assumption would be a' = 0. More precise estimates can be obtained 
from self-consistency arguments or from comparison with simulations. Both give 
a' /a £ 0.1. In the following we shall compare simulation data with predictions 
for the two extreme choices a' = and a' = 0.1a, mainly to show that the results 
depend very little on a'. 

For a' = our mean field treatment would be completely equivalent to percolation 
on a Bethe lattice with coordination number four, the sites of which are occupied 
with probability a (except for the root, which is occupied with probability g). This 
problem is exactly soluble and well discussed in the literature J7|. The extension 
to a' ^ is a straightforward application of the theory of branching processes ||. 

But first we have to compute g and a as functions of the controll parameter p. For 
this we use stationarity and the fact that each toppling decreases %i) by 4e, 
while each event (whether it actually triggers off an avalanche or not) increases 
it by one unit. Thus the average number of topplings per event is (s(p)) = l/4e. 
On the other hand, the average number of topplings during the first update of an 
avalanche is 4a, while it is multiplied by 3a + a' during each successive update. 
Taking into account that each event triggers an avalanche with probability g, we 
find 

{s{p)) = g[l + 4a£(3a + a')*] = g-—^ . (3) 

~^ 1 — 3a — a' 

Combining these two estimates, we find a somewhat complicated expression for g 
as a function of p. Instead of writing it down we just give inversely g and p as 
functions of a and a 1 : 

1 a — a' 

g= + 2a— - p = ag. (4) 

2 1 + a — a' 

This gives always p < 1/2, as we should expect since no stationary solution exists 
for p > 1/2. For p - l/2_, we find g = ^[1 - e ^fj+f^ ]. Fig.2 
shows the very reasonable agreement of this prediction with values obtained by 
simulations. We point out in particular that the data show that dg/dp is finite at 
p -»• 1/2, g w 0.6483 - 0.73e, as predicted by eq.(f§. 



Figure 2: Density g{p) against p. The continuous line is the mean field prediction of 
e Q-® with a' = 0, the dashed line is the prediction with a' = 0.1a, and the points show 
the numerical data obtained from simulations. 
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The quantity easiest to compute is the survival probability P t , defined as the 
probability that an event triggered an avalanche which lasts for > t time steps. We 
denote by Qt the probability that all sites are again stable at time t, provided that 
the event started with an unstable root at time 0. Then obviously P t = g(l — Qt)- 
Similarly, we call Qt the probability that an unstable site different from the root 
will not create any unstable offspring t time steps later. 

Furthermore, we denote by pk = — o) AL ~ k a k the probability that the first 

toppling triggers off k topplings at the next time step, and by pk the analog 
distribution for a later toppling. Their generating functions are given by 

4 

g(s) = Y J S k Pk = {l-a + saf (5) 

k=0 

and 

4 

g( s ) = s k p k = (1 - a + sa')(l -a + sa) 3 . (6) 
Then we have obviously 

Qt = T,Pk[Qt-i} k = g(Qt-i) (7) 

k 

and 

Pt = Q {l-Qt) = Qg{Qt-i). (8) 
In the vicinity of the critical point, Pt obeys the scaling law [[7] 

Pt « (9) 

with 5 = u t = 1. 



Figure 3: Avalanche survival probabilities Pt for different values of e. From left to right 
e = 0.1,0.0316,0.01,0.00316,0.001, and 0.000316. Solid lines are from mean field theory 
with a' = 0, dashed lines from mean field theory with a' /a = 0.1, and dots are from 
simulations. 

Predictions for P t are compared to simulation data in fig.^|. Notice that we used 
there the exact predictions from the recursion relation eq.(|7D which can be made 
numerically stable by some minor rearrangements. For small values of t we see 
perfect agreement, while there are deviations at very large t. They indicate that 
v t is somewhat larger than its mean field value, v t = 1.024 ± 0.008, while 5 = 
1.00 ± 0.01 is in exact agreement with the prediction. 

The other quantity we studied is the avalanche size distribution V n , defined as the 
probability that an event involves exactly n topplings. For n = and n = 1 we 
have Vo = 1 — f? and V\ = g(l — a) 4 . To calculate V n for n > 1 we use a theorem 
due to Dwass||. Consider a branching process where the distribution of offsprings 
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of a single individuum is given by pk, with generating function g(s). Its progeny 
is the total number of descendants, including itself. Then the probability that the 
total progeny of n individua is n is given by 

Pn\n = —Pn-no , U > n , (10) 

where is the k-th Taylor coefficient of [#(s)] n , i.e. [g(s)] n = p^ + PiS + 

p\ ' S l + 

Denoting again quantities referring to side branches by tilde's, we have in the 
present case 

4 4 k 

Vn = PkPn-l\k = 0J2 ——rPkPn-k-l , i 11 ) 

k=0 fc=0 71 1 

while p£ can be computed in several ways. The most straightforward is to use 
binomial expansions in [g(s)] ra . For large values of n we can also use the asymptotic 
behaviour. From the cumulant expansion of g(s) we see that p^ tends for n — > oo 
towards a Gaussian probability distribution with mean value n(k) and variance 
na, where (k) = 3a + a' and a = 3a(l — a) + a'(l — a') are the first two cumulants 
of p. This gives us 

no [(l-(fc))n-n ] 2 

Pn\n « ^==e W . (12) 

V 2Tin A a 



Figure 4: Integrated avalanche size distribution D n for the same values of e as in fig. 3. 
Solid lines are again predictions with a = 0, dashed lines are from predictions with 
a' = 0.1a, and points are from simulations. 

The integrated distribution D n = J2m= n P™, is shown in fig. 4. Again, the contin- 
uous and dashed lines give the predictions with a' = and a' = 0.1a, while the 
dots are results from simulations. This time the agreement is essentially perfect. 
This proves also that the simulation data satisfy the mean field scaling behaviour 



obtained from eq. (12 



D n = —V(ne 2 ) . (13) 

>n 



But again we found that the agreement with the detailed prediction of the Bethe 
lattice model is much better than that with the scaling form. 

Finally, we mention that we also performed some simulations at p = 1/2, on 
lattices of size L x L with open boundary conditions. We found that the average 
avalanche size (n) scales expected for any model of this type in which 

{J2i z i) is conserved 0, [10|] . 

In summary, we have shown numerically that a non-equilibrium model inspired 
by sandpile models shows a critical behaviour which is extremely close to mean 
field type. Indeed it seems that all exponents except one are identical to those for 
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mean field percolation. This is so in spite of the fact that the model involves only 
short range interactions and lives on a 2-d lattice. We have argued that this is 
due to the "confusion" brought about by the non-positive and non-linear effect of 
feedback loops. 

We should mention that similar behaviour was observed earlier in similar models 
fU . But apart from being more complicated, these models were still closer in spirit 
to the original sandpile model of [|1|]. In particular, in models were studied were 
the avalanche evolution was either stochastic or non- conservative, but not both 
together. It was the combination of both aspects which allowed us to draw a very 
close connection to percolation, and to compare in detail with percolation on a 
Bethe lattice. 

Finally, we should mention that our model is very similar - superficially seen - 
to the Ising model. Indeed, only very minor modifications were needed in the 
computer program given in H to switch from one model to the other. The main 
difference with the Ising model is that the present model does not satisfy detailed 
balance and is thus an inherently non-equilibrium system. We believe that the 
latter is necessary to observe mean-field type behaviour generated by the above 
mechanism. 
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